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ABSTRACT 

We present three-dimensional hydrodynamic simulations of the evolution of selfgrav- 
itating, thick accretion discs around hyperaccreting stellar-mass black holes. The 
black hole-torus systems are considered to be remnants of compact object mergers, 
in which case the disc is not fed by an external mass reservoir and the accretion is 
non-stationary. Our models take into account viscous dissipation, described by an a- 
law, a detailed equation of state for the disc gas, and an approximate treatment of 
general relativistic effects on the disc structure by using a pseudo-Newtonian poten- 
tial for the black hole including its possible rotation and spin-up during accretion. 
Magnetic fields are ignored. The neutrino emission of the hot disc is treated by a 
neutrino-trapping scheme, and the ^'t'-annihilation near the disc is evaluated in a 
post-processing step. Our simulations show that the neutrino emission and energy 
deposition by i/i^-annihilation increase sensitively with the disc mass, with the black 
hole spin in case of a disc in corotation, and in particular with the a- viscosity. We find 
that for sufficiently large a-viscosity j/i>-annihilation can be a viable energy source for 
gamma-ray bursts. 

Key words: accretion, accretion discs; black hole physics; gamma-rays: bursts; hydro- 
dynamics; neutrinos 
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1 INTRODUCTION 

In a series of papers (Ruffert, Janka & Schafer 1996; Ruffert 
et al. 1997; Ruffert & Janka 1999, 2001) it was shown that 
the neutrino emission associated with the dynamical phase 
of the merging or collision of two neutron stars (NS-f NS) is 
powerful, but too short to provide the energy for gamma- 
ray bursts (GRBs) by neutrino-antineutrino annihilation. 
Significant heating of the coalescing stars occurs only after 
they have plunged into each other. The neutrino luminosi- 
ties can then rise to several 10^"^ergs~^ (see also Rosswog 
& Liebendorfer 2003) and even exceed lO^^ergs"^ in case 
of the more violent (although probably very rare) collisions. 

After a few milliseconds, the compact massive rem- 
nant of the merger will most likely collapse to a black hole 
(BH) (see, e.g., Popham, Woosley & Fryer 1999, Di Mat- 
teo, Perna & Narayan 2002, OechsUn et al. 2004, Shibata & 
Uryu 2000). When this happens some matter can remain in 
a toroidal accretion disc around the BH and a funnel with 
low baryon density develops along the system axis. Hydrody- 
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namic simulations indicate that between some 0.01 Mq and 
a few 0.1 AIq of matter may have enough angular momen- 
tum to resist immediate absorption into the BH (see, e.g., 
Ruffert & Janka 1999, Shibata & Uryu 2000). A similar re- 
sult was obtained for BH-I-NS mergers (Janka et al. 1999, 
see also Lee 2001 and references therein). This matter is 
swallowed by the BH on the time scale of viscous transport 
of angular momentum, which is much longer than the dy- 
namical time scale. Different from the coUapsar case there 
is no external reservoir of stellar matter which feeds the 
accretion torus. Therefore the relevant time scale is set by 
viscosity-driven accretion rather than infall. The accretion is 
only approximately stationary, if global instabilities do not 
play a role (Ruffert et al. 1997, Lee & Ramirez-Ruiz 2002). 
Post-merger BH accretion might explain short GRBs, but is 
unlikely to account for long GRBs (e.g., Narayan, Piran & 
Kumar 2001). 

Post-merger accretion tori have maximum densities be- 
tween some 10^° gcm~^ and more than lO'^^ gcm~^ and ex- 
tend to outer radii of 10-20 times the Schwarzschild radius 
{Rb = 2GM/(?). The accretion proceeds with rates of frac- 
tions of a solar mass per second up to several solar masses 
per second. In case of such "hyperaccreting" black hole sys- 
tems (Popham et al. 1999) photons are strongly coupled to 
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the torus plasma and therefore are inefficient in transport- 
ing away energy. Neutrinos, however, are abundantly cre- 
ated by weak interactions in the very dense and hot tori 
and a fair fraction of the gravitational binding energy of the 
accreted matter can be radiated away by them ("neutrino- 
dominated accretion flow", NDAF, Popham et al. 1999, Ruf- 
fert et al. 1997, Ruffert & Janka 1999, Narayan et al. 2001, 
Kohri & Mineshige 2002, Di Matteo et al. 2002). 

In this Letter we present, for the first time, three- 
dimensional (3D) simulations of time-dependent hyperac- 
cretion from thick tori on a stellar-mass BH including the 
physical effects of on a stellar-mass BH including the physi- 
cal effects of BH rotation (Artemova, Bjornsson & Novikov 
1996), viscosity (a-prescription following Shakura & Sun- 
yaev 1973), a realistic finite-te mperat ure equation of state 
according to Lattimer & Swesty l|1991rl . and energy loss and 
change of lepton number by neutrino emission. The annihi- 
lation of emitted neutrinos and antineutrinos is investigated 
for its potential to provide the energy of relativistic gamma- 
ray burst fireballs (see, e.g., Eichler et al. 1989; Narayan, 
Paczyhski & Piran 1992; Woosley 1993). Our 3D modeling 
of the time-dependent (non-stationary) torus evolution abol- 
ishes approximations of previous (semi-)analytic work (e.g., 
Popham et al. 1999, Di Matteo et al. 2002) or more radically 
simplified numerical modeling in 2D (Lee & Ramirez-Ruiz 
2002). 



2 NUMERICS AND INITIAL CONDITIONS 

We use an Eulerian, 3D hydrodynamics code based on the 
Piecewise Parabolic Method of Colella & Woodward (1984). 
The code employs three nested Cartesian grids in a com- 
putational volume of 500 km side length. Each coarser grid 
level has the same number of zones but twice the zone size of 
the level below. Mirror symmetry relative to the equatorial 
plane is assumed. For an equatorial length and width of the 
computational volume of 500 km and a vertical extension of 
125 km, the smallest zones have a side length of 1.95 km. 

The gravitational potential of the black hole is described 
by the Artemova-Bjornsson-Novikov potential (Artemova et 
al. 1996), which reduces to the Paczynski-Wiita potential 
(Paczyhski & Wiita 1980) when the BH spin parameter 
a = Jc/{GM^ii) is zero. This potential mimics some gen- 
eral relativistic effects like the existence of an innermost 
stable circular orbit (ISCO) and the variation of the ISCO 
and event horizon with the value of a. The BH is repre- 
sented as a gravitating "vacuum sphere" in our calculations. 
The corresponding inner boundary radius is defined by the 
arithmetic average of event horizon and ISCO. Selfgravity 
of the torus is taken into account by a Newtonian descrip- 
tion. We us e the a-prescription proposed by Shakura & Sun- 
yaev lll973h to model the physical effects of disc viscosity 
(in addition to the ever present numerical viscosity), includ- 
ing the terms for viscous transport of angular momentum 
and dissipation of energy (a detailed description of the tech- 
nical implementation will be given elsewhere; Setiawan et 
al., in preparation). With the numerical resolution which 
we could afford in the present set of simulations, the dis- 
sipative effects of numerical viscosity should be somewhat 
larger than estimated for our previous work (see Ruffert & 
Janka 1999 and Janka et al. 1999) and might correspond to 



an Q- viscosity around 0.01. The latter value therefore sets a 
lower bound to the range of a's where we expect to observe 
influence by the inclusion of viscosity terms in the simu- 
lations discussed here. Gravitational-wave emission and its 
corresponding back-reaction (Blanchet, Damour & Schafer 
1990) are taken into account. A neutrino leakage scheme 
(Ruffert et al. 1996) is used to calculate the energy loss and 
lepton number change by neutrino emission from the disc 
gas. 

The initial configurations in our simulations consist of a 
BH surrounded by a toroidal accretion disc. These systems 
are considered as remnants of BH-fNS mergers and the sys- 
tem parameters were chosen appropriately. For this purpose 
we were guided by data obtained in some of the model runs 
described in Janka et al. (1999). The simulations start with 
a BH that has a mass of 4.017 Mq and an initial rotational 
parameter of or 0.6, respectively. Our "reference accretion 
torus" has a mass of 0.0478 Mq. This mass was changed 
to 0.0120 Mq or 0.19121/© by simply scaling the density 
distribution with appropriate factors. The initial tempera- 
ture is around 2MeV in the inner part of the torus (out 
to an equatorial distance from the system axis of ~ 70 km) 
and still above 1 MeV in the more dilute outer parts up to 
~ 130 km. The different combinations of model parameters 
of our simulations can be found in Tabled Varing the system 
parameters by scaling the initial torus density and setting 
the BH spin parameter might appear ad hoc. One should no- 
tice, however, that a BH-I-NS merger is an extremely violent 
event which should not be expected to produce a BH-torus 
configuration in perfect rotational equilibrium. Moreover, 
the torus properties will develop and adjust to the effects 
of viscous shear. This initial phase of relaxation after the 
start of our simulations does not last very long and most 
transients have died out already after ~ 20 ms. 



3 RESULTS 

The models listed in Table provide information about the 
influence of the torus mass and viscosity, and of the BH 
rotation on the system evolution. 

The compact NDAFs considered here advect most of 
the disk mass into the BH, and convective currents do not 
play an important role for mass and angular momentum 
transport (cf. Narayan et al. 2001). Angular momentum is 
carried outward mainly by shear effects. Because of the large 
ratio of BH to disk mass, the spin parameter a of the BH 
shows only minor variation in response to the angular mo- 
mentum that is associated with the mass swallowed by the 
BH. The value at the end of our computed evolution, af, 
has increased for a corotating disc, otherwise has decreased 
relative to ai (see Table 0. 

A comparison of Models rOO-64 and al3-64, which differ 
only by the presence of physical viscosity in the latter model, 
confirms that a = 0.01 is close to the threshold where the ef- 
fects of viscous terms become noticable. Except for the tem- 
perature, which is increased by shear heating, both models 
have very similar characteristic properties at the end of the 
simulations. Model al3-64 with the higher temperature has 
a correspondingly larger neutrino emission and is somewhat 
more inflated due to thermal pressure. The mass accretion 
rate is therefore marginally lower and the remaining torus 
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Table 1. Parameters and some results of the torus evolution models. In all cases the initial BH mass was Mgjj = 4.017 Mq, the number 
of zones of each level of the three nested grids was 64, and the simulation was performed over a time interval of Atcai = 40 ms. 
and ai are the gas mass and BH spin parameter at the beginning of the simulations. The direction of the BH rotation is indicated by 
"pro" and "ret", for prograde and retrograde relative to the disc, respectively. The physical viscosity parameter is a. All the following 
quantities are given at the end of the simulations: BH angular momentum parameter af, torus mass Mj, BH mass accretion rate Mj, 
typical accretion time scale of the torus, tacc = M^/M^, average torus density (p^) and average torus temperature (Tj) for the bulk of 
the torus gas, total neutrino luminosity, L^, and integral rate of energy deposition by i^P-annihilation around the accretion torus, -B^j^. 
The last column gives the total energy deposition by i/P-annihilation, Ej^i^ , in the time interval At^ai . 
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Figure 1. Black hole mass accretion rate, Md (left), total neutrino luminosity, L,, (middle), and integral rate of energy deposition by 
i/P-annihilation, Ej^i^ (right), as functions of time for different models. Note that only a small fraction (~ 1%) of this energy is released 
in the low-density funnel above the poles of the BH. 



after 40 ms is insignificantly more massive, in contrast to the 
consequences of a much higher viscosity (see below). 

Concerning the simulated torus evolution general trends 
can be observed in dependence of the varied parameters (i.e., 
Md, Oi, and a). The numbers given in Table Q reveal the 
following behaviour: 

• With higher disk mass Md the mass accretion rate Md 
increases, but the accretion time scale, iacc = Md/Md, shows 
Httle change. This is true, however, only for low disc viscosity 
and nonrotating BHs (compare Models rOO-64, irl-64, ir4- 
64). In case of high torus viscosity and spinning BHs (a = 
0.1; Models arl-64 and ar2-64) a four times larger initial 
disk mass correlates with a 5.6 times longer lifetime of the 
torus. 

• A BH in corotation with the disc leads to a higher torus 
temperature, because the ISCO and event horizon shrink 
and the torus moves closer to the BH. Since thermal pressure 
puffs up the torus (as visible in Fig.|2l, its density does not 



rise in the same way. The mass accretion rate decreases and 
the accretion time scale grows. Therefore the torus mass is 
higher at a given time. Retrograde rotation of the BH has 
the opposite effect (compare Models rOO-64, ro2-64, ro5-64 
or al4-64 with arl-64). 

• Without physical viscosity, the tori remain relatively 
cool (near their initial temperature), independent of the 
torus mass. Larger viscosity (bigger a) increases the mass 
accretion rate and reduces the accretion time scale face- This 
holds for nonrotating BHs (compare Models rOO-64 and al4- 
64) and rotating BHs (Models ro2-64 and arl-64). With 
higher torus viscosity the temperature becomes higher and 
the neutrino emission correspondingly stronger. Torus mass 
and density are lower at the same time. The similarity of 
Models al4-64 and arl-64, and sizable differences relative to 
Model al3-64 suggest that increasing the viscosity from 0.01 
to 0.1 is more significant than BH rotation with a = 0.6 
instead of no rotation. 
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Figure 2. Efficiency of conversion of rest-mass energy of accreted 
matter to neutrino emission, = / {Mc^) (tfiick lines), and 
j/iz-anniiiilation efficiency, q^^ = Euj^/Lu (tliin lines), for different 
models as functions of time. 

The total neutrino luminosity of the torus, L^, (i.e., the 
sum of the luminosities of neutrinos and antineutrinos of 
all flavors) increases with the torus mass, with the BH spin 
in case of direct rotation (corotation with disc), and with 
the viscosity. This is a consequence of a higher torus tem- 
perature. Again the sensitivity is largest to a variation of 
the viscosity within the considered bandwidth of values. For 
each model phases of high neutrino luminosity correlate with 
large mass accretion rates (compare the left and middle pan- 
els of Fig. 0, indicating that the processes which lead to a 
high value of Aid (in particular during the first 20 ms of the 
computed evolution) also cause heating of the torus plasma. 
But this correlation is not strict. This can be seen in case 
of Model ir4-64 where after ~ 15 ms the mass accretion rate 
becomes nearly constant but shows a continuous trend 
of increase. 

Since the integral rate Ei^c of energy deposition by 
neutrino-antineutrino annihilation to electron-positron pairs 
(ui? — > e^e~) in the surroundings of the BH-torus system 
scales with the product of neutrino and antineutrino lumi- 
nosities, 

= CL.LJil%l±MlM) , (1) 

this rate follows closely the behavior of L^. This can be seen 
by a comparison of the middle panel of Fig. Q which shows 
Lv{t), with the right panel of this figure, which gives Evvit). 
The factors {e^) and (e^) in Eq. Q denote the mean energy 
and mean squared energy for and u. The factor C contains 
weak interaction coefficients and terms that depend on the 
geometry of the neutrino-emitting torus region. The latter 
does not change much during the accretion process. Note 
that in the considered problem only Ve and need to be 
taken into account for computing the annihilation rate. The 



contributions from muon and tau neutrinos and antineu- 
trinos can safely be ignored because their luminosities are 
lower than those of Ve and Ue by factors of several, and their 
annihilation cross section is also smaller by a factor of about 
5. The general decline of E^v that is visible after 15-20 ms 
in all cases except Model ir4-64, can be rather well fitted by 
a t~^^'^ power-law. 

Figure |5| displays the conversion efficiency of rest-mass 
energy to neutrinos, = / [Mc?), as function of time, 
and the corresponding efficiency for conversion of neutrino 
energy to e^e~-pairs by i/j/-annihilation, q^y = Evj^/Lu- 
Values of several per cent can be reached for both quan- 
tities in case of the models with viscosity a — 0.1, Models 
al4-64, arl-64, and ar2-64. 

From the results listed in Table and Figs. and |5| 
it is clear that only models with large disc viscosity pro- 
duce sufficiently high neutrino luminosities to allow for 
!^/?-annihilation as potential energy source of cosmological 
GRBs. This is apparent from Fig.|3where ray-tracing images 
of the neutrino emission of the non-viscous Models rOO-64, 
ro2-64, and ro5-64 are contrasted by the high-viscosity Mod- 
els al4-64, arl-64, and ar2-64. The latter models show high 
temperature and intense neutrino emission in a much more 
extended volume around the BH. Despite of the different 
impression left by the lower plots of Fig. |3] there is only one 
of the investigated models where the torus is optically thick 
to neutrinos during the computed evolution, namely Model 
ir4-64 (not displayed in Fig. EJ, which has quite a massive 
torus with the highest density of all models in Tabled This 
confirms that in most cases the treatment of neutrino effects 
by a trapping scheme is adequate and neutrino transport is 
not important. 

Model ar2-64 is the most extreme case with high vis- 
cosity, high torus mass, and direct rotation of a Kerr-BH. It 
combines the conditions which are most favorable for pro- 
ducing GRBs. It shows the highest torus temperatures and 
neutrino luminosities, yielding a total energy deposition rate 
by !/i^-annihilation of more than 10^^ ergs~^ at 40 ms when 
the simulations were stopped. At f = 12 ms the value of E^a 
is even around 10^^ ergs~^ (Figs.Qand^J- Most of this en- 
ergy, however, is deposited close to the equatorial plane in a 
region where the torus density is high and the rate of energy 
loss by neutrino emission is larger than the energy input by 
z/iz-annihilation (cf. the middle and right panels of Fig. IIJ. 
Therefore this energy is useless for powering relativistic out- 
flow. But a low-density funnel has developed along the sys- 
tem axis above the poles of the BH. Within this funnel the 
energy deposition by i/iz-annihilation exceeds the local ree- 
mission of neutrinos and accounts for an energy input rate 
to that region of order lO^^ergs"^. With a torus lifetime 
of several tenths of a second (Table 0, some 10*^ erg might 
thus be suited to power a pair of ultrarelativistic axial e'^e~- 
plasma jets. The involved energy is sufficiently large to make 
i/i/-annihilation a viable energy source for short GRBs. 



4 DISCUSSION 

We have performed three-dimensional hydrodynamic sim- 
ulations of hyperaccreting stellar mass BHs with self- 
gravitating, thick accretion tori, varying the torus mass, 
BH spin and a-viscosity of the torus gas. Our models show 
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Figure 3. Ray-tracing images of the neutrino emission of the zero-viscosity Models rOO-64, ro2-64, and ro5-64 (upper plots, from left 
to right) in comparison to Models al4-64, arl-64, and ar2-64 (lower plots, from left to right) which were computed with a viscosity value 
of a = 0.1. The colors (ranging from red, orange, yellow, green, blue to white) visualize the gas temperature (roughly 2, 4, 6, 8, 10, and 
12MeV or higher), the brightness is a measure of the intensity of the neutrino emission. The figures are snapshots that correspond to 
times (from top left to bottom right) 5.5, 5.8, 5.5, 14.2, 7.0, and 14.0 ms, respectively, after the start of the simulations. A region with a 
diameter of about 500 km around the BH at the centre is displayed. The dark spheres have a radius given by the arithmetic mean of the 
event horizon and ISCO of the BH with values of, from top left to bottom right, 24, 17, 29, 24, 17, and 17 km, respectively. The white 
dots symbolize background "stars" which are obscured by a neutrino-opaque accretion torus. The zero-viscosity tori are sufficiently hot 
and dense only very close to the BH (in particular in case of direct BH rotation. Model r02-64), while viscosity leads to heating and high 
neutrino emission in a much more extended volume. 
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Figure 4. Neutrino emission and I'P-annihilation in Model ar2-64 at time 12.18 ms. The left panel displays the surface emissivity of 
Ue projected onto the equatorial plane. The contours are labeled with the logarithm of the energy loss rate in ergcm^^s"'^. They are 
spaced in steps of half a unit. The dotted lines indicate the region where the optical depth for Ue exceeds unity; the torus is essentially 
transparent to neutrinos. In the middle panel the energy loss rate per unit volume (in ergcm~''s~^) in neutrinos and antineutrinos 
of all flavors is plotted in the x-z- and y-z-planes perpendicular to the equatorial plane. The contours represent the logarithm of the 
rate and are spaced in steps of half a unit. The right panel shows a map of the local energy deposition rate (in ergcm~^s~^) by uu 
annihilation in the surroundings of the BH which is indicated by the white octagon at the centre. The contours are logarithmically 
spaced with the darker grey shading meaning higher energy deposition rate. The integral value of the deposition rate at the displayed 
time is 8.5 X 10^^ ergs~^. The dotted contour lines represent levels of constant values of the azimuthally averaged mass density (spaced 
logarithmically with intervals of 0.5 dex) with the bold dotted line corresponding to p = 10^" gcm~'^. 
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that tori which are heated by viscous dissipation of energy 
get inflated by thermal pressure and are therefore unlikely 
to become optically thick to neutrinos. The neutrino lu- 
minosities stay well below the Eddington limit L^^Edd = 
4-kGMbhc/k^ ~ 2 X 10^^ (MBH/4Mo)ergs"^ Here ~ 
10-"(e^)/(20MeV)^ is the mean neutrino opacity (e.g., 
Janka 2001). The latter is mostly determined by the con- 
tribution of S-e which dominate the energy loss of the torus. 
These findings are in contrast to conclusions drawn from 
more approximative modeling approaches which treated the 
vertical structure of the torus by height- averaging (Di Mat- 
teo et al. 2002). 

Our results also suggest that the integral rate of the 
energy deposition by i^i/-annihilation drops like t~^^^ dur- 
ing the long-time, slow decay of the accretion rate that fol- 
lows an initial, transient relaxation phase of 10-20 ms with 
very high mass accretion rates and neutrino emission. This 
temporal decrease is less steep than obtained in previous 
simulations with azimuthal symmetry which used a simple 
ideal gas equation state and made the assumption that all 
the dissipated energy is radiated away in neutrinos (Lee & 
Ramirez- Ruiz 2002). 

We have shown that :/i?-annihilation in the low-density 
funnel above the poles of the BH can deposit energy at 
a rate of ~ lO^^ergs"^, accounting for a total energy re- 
lease of some 10*^ erg in case of sufficiently large torus mass 
(> 0.1 M0), high disk viscosity (a ~ 0.1) and BH rota- 
tion with spin parameter a ~ 0.6. These conditions are ex- 
pected to be generically produced by BH-I-NS mergers and 
most likely also by NS-I-NS mergers (Janka et al. 1999). 
The main effect of direct BH rotation is an increase of the 
lifetime of the torus. The energy release by z/iz-annihilation 
satisfies the energy requirements of short GRBs in case of 
a moderate amount of collimation of the ultrarelativistic 
bipolar outfiows into a fraction / — 2A^l/4n of a few per 
cent of the sky (e.g., Rosswog & Ramirez-Ruiz 2003). Our 
results for the post-merging accretion of a remnant BH- 
torus system are therefore more optimistic than the esti- 
mates based on NS-I-NS merger simulations by Rosswog & 
Ramirez- Ruiz (2002). Future simulations, however, will have 
to show whether the deposited thermal energy can be effi- 
ciently converted to axial outffows from BH-torus systems 
that are sufficiently luminous and coUimated to account for 
short GRBs (Aloy, Janka & Miiller, in preparation). 
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